home *** CD-ROM | disk | FTP | other *** search
/ ADA Programming Guide / ADA Programming Guide.iso / ada_gwu / adaed / math / primfunc.adb < prev    next >
Encoding:
Text File  |  1996-01-30  |  9.4 KB  |  335 lines

  1. --  -----------------------------------------------------------------------
  2. --  Title:       primitive_functions
  3. --  Last Mod:    Thu Apr 19 11:25:43 1990
  4. --  Author:      Vincent Broman <broman@nosc.mil>
  5. --   Copyright 1990 Vincent Broman
  6. --     Permission granted to copy, modify, or compile this software for
  7. --   one's own use, provided that this copyright notice is preserved intact.
  8. --     Permission granted to distribute compiled binary copies of this
  9. --   software which are linked in with some other application.
  10. --     Permission granted to distribute other copies of this software,
  11. --   provided that (1) any copy which is not source code, i.e. not in the
  12. --   form in which the software is usually maintained, must be accompanied
  13. --   by a copy of the source code from which it was compiled, and (2) the
  14. --   one distributing it must refrain from imposing on the recipient
  15. --   further restrictions on the distribution of this software.
  16. --   
  17. --  Visibility:  
  18. --  Description: 
  19. --               functions on floating point types whose implementation
  20. --               depends on access to the mantissa/exponent representation
  21. --               of the floating point number.  this includes
  22. --               integer/fraction operations.
  23. --               
  24. --               This is a portable, slow implementation.
  25. --               
  26. --  Exceptions:  numeric_error upon overflow of function scale,
  27. --               constraint_error only if type Float is constrained.
  28. --  -----------------------------------------------------------------------
  29. package body primitive_functions is
  30. -- 
  31.    
  32.    
  33.    subtype expbits is integer range 0 .. 6; -- 2**(2**7) might overflow
  34.  
  35.    Float_exp: array( expbits) of integer :=
  36.       (1, 2, 4, 8, 16, 32, 64);
  37.    Float_power: array( expbits) of Float :=
  38.       (2.0 **  1, 2.0 **  2, 2.0 **  4, 2.0 ** 8,
  39.        2.0 ** 16, 2.0 ** 32, 2.0 ** 64);
  40.    Float_neg_power: array( expbits) of Float :=
  41.       (0.5 **  1, 0.5 **  2, 0.5 **  4, 0.5 ** 8,
  42.        0.5 ** 16, 0.5 ** 32, 0.5 ** 64);
  43.  
  44.   -- double stuff deleted here 
  45.  
  46.    function exponent( x: Float) return integer is
  47. --
  48. -- return the exponent k such that 1/2 <= x/(2**k) < 1,
  49. -- or zero for x = 0.0 .
  50. -- 
  51.    begin
  52.       if x = 0.0 then
  53.      return 0;
  54.       end if;
  55.       -- nonzero x
  56.       -- essentially, just multiply by powers of two to get in range
  57.       
  58.      declare
  59.         ax: Float := abs(  x);
  60.         exp: integer := 0;
  61.         -- ax*2.0**exp is invariant
  62.      begin
  63.         while ax >= Float_power( expbits'last) loop
  64.            ax := ax * Float_neg_power( expbits'last);
  65.            exp := exp + Float_exp( expbits'last);
  66.         end loop;
  67.         -- ax < 2**64
  68.         for n in reverse expbits'first .. expbits'last - 1 loop
  69.            if ax >= Float_power( n) then
  70.           ax := ax * Float_neg_power( n);
  71.           exp := exp + Float_exp( n);
  72.            end if;
  73.            -- ax < Float_power( n)
  74.         end loop;
  75.         -- ax < 2
  76.         while ax < Float_neg_power( expbits'last) loop
  77.            ax := ax * Float_power( expbits'last);
  78.            exp := exp - Float_exp( expbits'last);
  79.         end loop;
  80.         -- 2**-64 <= ax < 2
  81.         for n in reverse expbits'first .. expbits'last - 1 loop
  82.            if ax < Float_neg_power( n) then
  83.           ax := ax * Float_power( n);
  84.           exp := exp - Float_exp( n);
  85.            end if;
  86.            -- Float_neg_power( n) <= ax < 2
  87.         end loop;
  88.         -- 1/2 <= ax < 2
  89.         if ax < 1.0 then
  90.            return exp;
  91.         else
  92.            return exp + 1;
  93.         end if;
  94.      end;
  95.       
  96.    end exponent;
  97.  
  98.  
  99.    function mantissa (x: Float) return Float is
  100. --
  101. -- return scale( x, - exponent( x)) if x is nonzero,  0.0 otherwise.
  102. --
  103.    begin
  104.       if x = 0.0 then
  105.      return 0.0;
  106.       end if;
  107.       -- nonzero x
  108.       -- essentially, just multiply by powers of two to get in range
  109.       
  110.      declare
  111.         ax: Float := abs( x);
  112.      begin
  113.         while ax >= Float_power( expbits'last) loop
  114.            ax := ax * Float_neg_power( expbits'last);
  115.         end loop;
  116.         for n in reverse expbits'first .. expbits'last - 1 loop
  117.            if ax >= Float_power( n) then
  118.           ax := ax * Float_neg_power( n);
  119.            end if;
  120.         end loop;
  121.         while ax < Float_neg_power( expbits'last) loop
  122.            ax := ax * Float_power( expbits'last);
  123.         end loop;
  124.         for n in reverse expbits'first .. expbits'last - 1 loop
  125.            if ax < Float_neg_power( n) then
  126.           ax := ax * Float_power( n);
  127.            end if;
  128.         end loop;
  129.         if ax >= 1.0 then
  130.            ax := ax * 0.5;          -- 1 <= ax < 2
  131.         end if;
  132.         if x > 0.0 then
  133.            return Float( ax);
  134.         else
  135.            return Float( - ax);
  136.         end if;
  137.      end;
  138.       
  139.    end mantissa;
  140.       
  141.  
  142.    function scale (x: Float;
  143.            k: integer) return Float is
  144. --
  145. -- return x * 2**k quickly, or quietly underflow to zero,
  146. -- or raise an exception on overflow.
  147. -- 
  148.    begin
  149.       if x = 0.0 or k = 0 then
  150.      return x;
  151.       end if;
  152.       -- nonzero x and k.
  153.       -- essentially, just multiply repeatedly by 2**(+-2**n).
  154.       -- if is_
  155.      declare
  156.         y: Float  :=  x;
  157.         exp: integer := k;
  158.         -- y*2.0**exp is invariant
  159.      begin
  160.         while exp <= - Float_exp( expbits'last) loop
  161.            y := y * Float_neg_power( expbits'last);
  162.            exp := exp + Float_exp( expbits'last);
  163.         end loop;
  164.         for n in reverse expbits'first .. expbits'last - 1 loop
  165.            if exp <= - Float_exp( n) then
  166.           y := y * Float_neg_power( n);
  167.           exp := exp + Float_exp( n);
  168.            end if;
  169.         end loop;
  170.         -- exp >= 0
  171.         while exp >= Float_exp( expbits'last) loop
  172.            y := y * Float_power( expbits'last);
  173.            exp := exp - Float_exp( expbits'last);
  174.         end loop;
  175.         for n in reverse expbits'first .. expbits'last - 1 loop
  176.            if exp >= Float_exp( n) then
  177.           y := y * Float_power( n);
  178.           exp := exp - Float_exp( n);
  179.            end if;
  180.         end loop;
  181.         -- exp = 0
  182.         return Float( y);
  183.      end;
  184.       
  185.    end scale;
  186.  
  187. --
  188. -- truncation requires making all fractional bits zero in a signed
  189. -- magnitude representation of floating point numbers.
  190. -- for large numbers, this has no effect, as all large floats are integers.
  191. -- for abs( x) < 1.0 this produces zero.
  192. --
  193.  
  194.    function truncate (x: Float) return Float is
  195. --
  196. -- truncate x to the nearest integer value with absolute value
  197. -- not exceeding abs( x).  No conversion to an integer type
  198. -- is expected, so truncate cannot overflow for large arguments.
  199. --
  200. -- 
  201.       large: Float  := 1073741824.0;
  202.       type long is range - 1073741824 .. 1073741824;
  203.       -- 2**30 is longer than any single-precision mantissa
  204.       rd: Float;
  205.  
  206.    begin
  207.       if abs( x) >= large then
  208.      return x;
  209.       else
  210.      rd := Float ( long( x));
  211.      if x >= 0.0 then
  212.         if rd <= x then
  213.            return rd;
  214.         else
  215.            return rd - 1.0;
  216.         end if;
  217.      else
  218.         if rd >= x then
  219.            return rd;
  220.         else
  221.            return rd + 1.0;
  222.         end if;
  223.      end if;
  224.       end if;
  225.    end truncate;
  226.  
  227.  
  228.    function shorten (x: Float;
  229.              k: integer) return Float is
  230. -- 
  231. -- set all but the k most significant bits in the mantissa of x to zero,
  232. -- i.e. reduce the precision to k bits, truncating, not rounding.
  233. -- shorten( x, k) = 0.0 if k < 1 and
  234. -- shorten( x, k) = x   if k >= Float'machine_mantissa.
  235. -- 
  236.       kex: constant integer := k - exponent( x);
  237.  
  238.    begin                                -- shorten
  239.       -- the tests on k are needed only to avoid under/overflow
  240.       if x = 0.0 or k < 1 then
  241.      return 0.0;
  242.       elsif k >= Float'machine_mantissa then
  243.      return x;
  244.       else
  245.      return scale( truncate( scale( x, kex)), - kex);
  246.       end if;
  247.    end shorten;
  248.  
  249.  
  250.    function odd (x: Float) return boolean is
  251. --
  252. -- predicate indicates whether or not truncate( x) is an odd integer.
  253. --
  254.      x2: constant Float := truncate(x) / 2.0;
  255.  
  256.    begin
  257.       return x2 /= truncate( x2);
  258.    end odd;
  259.  
  260. --
  261. -- all the other integer/fraction functions are based on truncate.
  262. --
  263.  
  264.    function floor (x: Float) return Float is
  265. --
  266. -- return as a Float the greatest integer value <= x.
  267. --
  268.         sx: constant Float  := x; 
  269.         tx: constant Float := truncate( sx); -- no constraint_error
  270.      begin
  271.         if sx >= 0.0 or sx = tx then
  272.            return Float( tx);
  273.         else
  274.            return Float( tx - 1.0); -- exactly
  275.         end if;
  276.      
  277.       
  278.    end floor;
  279.  
  280.  
  281.    function ceiling (x: Float) return Float is
  282. --
  283. -- return as a Float the least integer value >= x.
  284. --
  285.         sx: constant Float :=  x;
  286.         tx: constant Float := truncate( sx); -- no constraint_error
  287.      begin
  288.         if sx <= 0.0 or sx = tx then
  289.            return Float( tx);
  290.         else
  291.            return Float( tx + 1.0); -- exactly
  292.         end if;
  293.      
  294.       
  295.    end ceiling;
  296.  
  297.  
  298.    function round (x: Float) return Float is
  299. --
  300. -- return as a Float the integer value nearest x.
  301. -- in case of a tie, prefer the even value.
  302. --
  303.    
  304.       
  305.         sx: Float  :=  x;
  306.         tx: Float := truncate( sx);
  307.         diff: Float := sx - tx;     -- computed exactly
  308.      begin
  309.         if abs( diff) < 0.5 then
  310.            return Float( tx);
  311.         elsif diff > 0.5 then
  312.            return Float( tx + 1.0);
  313.         elsif diff < -0.5 then
  314.            return Float( tx - 1.0);
  315.         elsif diff = 0.5 then
  316.            if odd( tx) then
  317.           return Float( tx + 1.0);
  318.            else
  319.           return Float( tx);
  320.            end if;
  321.         else                        -- diff = -0.5
  322.            if odd( tx) then
  323.           return Float( tx - 1.0);
  324.            else
  325.           return Float( tx);
  326.            end if;
  327.         end if;
  328.      
  329.       
  330.    end round;
  331.  
  332. end primitive_functions;
  333.  
  334. -- $Header: p_primitive_functions_b.ada,v 3.13 90/04/19 12:34:22 broman Rel $
  335.